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The  single-particle  model  is  a  useful  mathematical  representation  for  state-of-charge  observation,  param¬ 
eter  identification  and  control  of  lithium-ion  batteries.  This  model  is  a  simplified  electrochemical 
formulation  where  ionic  intercalation,  described  as  a  diffusive  process,  is  considered  as  the  dominant 
dynamics.  In  the  search  for  a  more  efficient  numerical  solution  of  the  involved  partial  differential  equa¬ 
tions,  many  approximations  have  been  reported.  However,  most  of  them  are  valid  just  under  restricted 
operating  conditions.  In  this  paper,  spatial  semidiscretization  is  reintroduced  as  the  precision  of  approx¬ 
imations  could  be  arbitrarily  chosen  with  this  approach.  Three  discretization  methods  are  applied  in  a 
classical  fashion,  evaluated  and  compared  in  both  time  and  frequency  domains:  finite  difference,  finite 
element  and  differential  quadrature.  In  addition,  two  commonly  used  low  order  approximations  are 
tested  against  semidiscretization  approximations.  The  best  results  are  obtained  with  the  differential 
quadrature  method  in  its  polynomial  version.  Two  model  truncation  criteria  are  also  explored,  one  is 
based  on  bandwidth  selection  and  the  other  on  residue  analysis,  where  the  first  resulted  to  be  more 
conservative.  Finally,  simulations  of  representative  reduced  order  approximations  of  the  single-particle 
model  are  compared  against  experimental  data  found  in  the  literature. 

©  201 1  Elsevier  B.V.  All  rights  reserved. 


1.  Introduction 

Due  to  a  natural  minimization  trend,  portable  applications  from 
cellular  phones  to  hybrid  electric  vehicles  require  more  efficient 
power  sources  such  that  energy  could  be  stored  in  smaller  and  more 
lightweight  devices.  On  the  other  hand,  renewable  resources  like 
wind  and  solar  energies  are  being  used  to  complement  or  substitute 
traditional  sources  of  electricity  generation,  helping  to  reduce  the 
production  of  pollutants  and  the  shortage  of  fossil  fuels.  However, 
as  alternative  energy  sources  depend  on  weather,  electricity  gener¬ 
ated  in  excess  should  be  stored  and  then  released  when  demand  is 
greater  than  production.  Thereby,  energy  storage  becomes  crucial, 
making  possible  the  existence  of  portable  appliances  and  trans¬ 
forming  renewable  sources  into  reliable  supplies. 

Reviews  about  storage  devices  currently  used  in  power  systems 
could  be  found  in  [1,2].  Electrochemical  batteries  have  represented 
the  most  popular  alternative  for  more  than  one  hundred  years. 
Historical  outlines  about  development  of  battery  technologies, 
comparisons  of  materials  and  chemistries  as  well  as  specifications 
can  be  found  in  [3-5].  Batteries  are  in  general  low-power  and  high- 
energy  storage  devices;  they  can  sustain  plain  charge  or  discharge 
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conditions  for  long  periods  but  are  not  well  suited  to  manage  power 
peaks.  Dominant  families  of  secondary  (reversible  or  rechargeable) 
batteries  are  the  lead-acid,  nickel-based  and  lithium-ion  (Li-Ion) 
systems.  This  last  family  represents  the  most  promising  technology, 
displacing  nickel-based  systems  in  the  fields  of  electronic  gadgets 
and  medical  implants,  and  is  expected  to  be  present  in  coming 
hybrid  electric  vehicles  ([6-8]). 

Li-Ion  batteries  posses  the  best  performance  features.  For 
lead-acid  batteries  specific  energy  and  power  are  30Whkg-1 
and  180  W  kg-1,  whereas  these  specifications  are  80-150  Wh  kg-1 
and  500-2000  W  kg-1  for  Li-Ion,  and  100-150  Wh  kg-1  and 
50-250  W  kg-1  for  Li-Ion  polymer  batteries,  respectively.  Li-Ion 
systems  do  not  present  the  memory  effect,  typical  of  the  nickel- 
based  family,  have  energy  efficiencies  of  90-100%,  self  discharge 
rates  as  low  as  5%  per  month,  can  reach  more  than  1500  full 
charge-discharge  cycles  and  have  an  open  circuit  potential  around 
3.7  V  [2].  Although  Li-Ion  batteries  have  low  environmental  impact, 
many  challenges  are  being  faced  in  regard  to  materials  and  manu¬ 
facturing  processes  to  make  this  technology  sustainable.  Research 
about  organic  compounds,  nanomaterials  and  lithium  alternatives 
like  aluminum  and  magnesium  are  being  strongly  impelled  [9]. 

Mathematical  models  of  batteries  are  combinations  of  static 
and/or  dynamic  relationships  that  set  a  correspondence  between 
current  and  voltage  at  battery  terminals  (impedance),  the  only  two 
measurable  variables.  Modeling  is  often  intended  to  on-line  esti¬ 
mation  of  non-measurable  variables  as  the  state-of-charge  and 
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the  state-of-health  (aging),  as  well  as  impedance  identification. 
Other  important  objective  is  off-line  battery  simulation  for  design. 
Simulations  make  it  possible  the  improvement  of  materials  and 
optimization  without  actual  manufacturing,  reducing  development 
time  and  costs  [3]. 

Li-Ion  battery  models  are  in  general  of  two  types:  equivalent  cir¬ 
cuits  and  electrochemical  models.  The  first  are  very  popular,  being 
intuitive  for  people  familiar  with  analysis  of  electric  networks. 
There  is  a  wide  variety  of  them  and  are  mainly  used  for  impedance 
identification.  Simple  circuits,  as  the  embedded  in  the  Advanced 
Vehicle  Simulator  (ADVISOR)  of  the  US  National  Renewable  Energy 
Laboratory  [10],  or  those  reported  in  [11,12],  are  extensively  used 
in  practice,  but  more  elaborated  versions  as  [13-15]  are  also 
available.  For  these  models,  constitutive  relationships  of  elements 
are  complicated  heuristic  functions  of  state-of-charge  (SOC)  and 
temperature,  and  their  use  for  analysis  results  difficult.  SOC  is 
the  fundamental  variable  for  battery  control.  Nevertheless,  when 
equivalent  circuits  are  used,  SOC  is  estimated  through  numerical 
integration  of  the  applied  current  (Coulomb  counting),  a  not  very 
accurate  even  robust  method  which  drives  to  conservative  power 
and  energy  management  algorithms. 

Electrochemical  models  were  originally  thought  for  battery 
design  and  optimization  as  they  offer  a  better  illustration  about  cell 
functioning  and  limitations.  The  formulation  is  based  on  physical 
and  chemical  principles,  resulting  in  sets  of  coupled  partial  dif¬ 
ferential  and  algebraic  equations.  This  approach  allows  accurate 
calculations  of  the  spatial  and  temporal  evolution  of  internal  vari¬ 
ables.  Basic  electrochemical  models  stated  for  charge-discharge 
prediction  of  Li-Ion  cells  with  two  porous  electrodes,  the  systems 
currently  used  for  power  applications,  can  be  found  in  [16-19], 
whereas  a  more  formal  reformulation  is  shown  in  [20].  Addition  of 
further  phenomena  to  the  basic  models  like  temperature  changes, 
electric  double-layer  capacitance  effect,  battery  aging  and  porosity 
variation  are  explored  in  [21-24],  for  example. 

Obtaining  rigorous  numerical  solutions  of  electrochemical  mod¬ 
els  is  time  consuming  and  requires  substantial  computational 
resources.  Equations  are  strongly  coupled  and  this  worsens  as 
supplementary  phenomena  are  appended.  Even  though,  because 
of  their  physical  meaning,  high  accuracy  and  more  realistic  SOC 
predictions,  many  efforts  are  being  carried  out  to  find  appropri¬ 
ate  simplifications  and  reduced  order  approximations  that  allow 
electrochemical  models  to  be  useful  for  on-line  estimation  and 
identification;  see  for  example  [25-33].  Depending  on  the  appli¬ 
cation  requirements,  the  first  step  is  to  define  which  dynamics  are 
of  interest,  for  example,  concentration  of  lithium  ions  in  the  solid 
state,  concentration  of  electrolyte  in  the  electrolytic  solution,  the 
effect  of  double  layer  capacitance,  cell  aging,  temperature  changes 
or  any  combination.  Then,  the  other  dynamics  could  be  either 
neglected  or  supposed  static,  driving  to  simpler  models.  As  most 
of  the  phenomena  involved  in  the  operation  of  a  lithium-ion  cell 
are  described  by  distributed  parameter  equations,  simplified  mod¬ 
els  can  contain  infinite  order  equations  that  should  be  efficiently 
numerically  solved  for  on-line  implementation. 

In  this  paper,  the  single-particle  model  is  studied  as  it  con¬ 
stitutes  a  useful  mathematical  representation  for  state-of-charge 
observation,  parameter  identification  and  control  of  lithium-ion 
batteries.  This  model  is  a  simplified  electrochemical  formulation 
where  ionic  intercalation,  described  as  a  diffusive  process,  is  con¬ 
sidered  as  the  dominant  dynamics.  Such  phenomenon  is  directly 
related  to  the  SOC  of  the  cell  and  voltage  variations  at  the  cell  ter¬ 
minals.  The  other  dynamics,  referred  above,  are  neglected  because 
they  have  a  weaker  influence  on  the  relationship  between  input 
current  and  output  voltage.  In  particular,  changes  in  temperature 
are  ignored  as  the  focus  in  this  paper  is  on  studying  the  effect  of 
continuous  charging-discharging  cycles  when  a  steady  sate  tem¬ 
perature  has  been  reached. 


The  state  equation  of  the  single-particle  electrochemical  model 
is  studied  with  the  goal  of  finding  reduced  order  approximations 
that  represent  adequately  the  battery  main  dynamics  for  some 
given  operation  conditions.  Spatial  semidiscretization  approach  is 
reintroduced  for  order  reduction  and  three  discretization  meth¬ 
ods,  combined  with  four  distributions  of  discretization  points, 
are  applied  in  a  classical  fashion:  finite  difference,  finite  element 
and  differential  quadrature.  Because  of  their  recurrent  appear¬ 
ance  in  the  literature,  dynamical  approximations  based  on  second 
and  fourth  order  polynomials  are  also  evaluated  and  compared 
against  semidiscretization  approximations.  Additionally,  two  order 
selection  or  truncation  criteria,  suggested  but  not  explored  in  the 
literature,  are  discussed  and  contrasted.  The  first,  more  conserva¬ 
tive,  consists  on  selecting  the  bandwidth  of  the  approximation  such 
that  it  covers  at  least  the  major  part  of  the  energy  spectrum  of  some 
proposed  test  signal.  The  second  criterion  is  more  relaxed  and  is 
based  on  a  residue  analysis  of  the  convolution  between  the  state 
equation  and  the  test  signal,  driving  to  more  practical  results. 

Approximations  are  first  evaluated  in  the  time  domain  to  have 
a  preliminary  surmise  about  which  is  more  accurate  in  comparison 
to  a  base  solution.  In  this  case,  the  test  signal  consists  on  a  rectangu¬ 
lar  pulse  followed  by  a  relaxation  time,  whose  duration  is  such  that 
allows  to  observe  the  whole  transient  response  under  forced  as  well 
as  under  free  excitation.  Then,  the  frequency  response  of  approx¬ 
imations  is  compared  against  that  of  an  analytic  transfer  function 
related  to  the  state  equation  in  order  to  confirm  the  results  obtained 
in  the  time  domain  tests.  Finally,  representative  cases  of  all  con¬ 
sidered  approximations  are  applied  to  the  single-particle  model 
and  their  response  to  the  test  signal  specified  in  the  FreedomCAR 
manual  [34]  is  compared  against  experimental  data  extracted  from 
[35].  Results  indicate  that  differential  quadrature  approximations 
are  the  best  and  that  low  order  polynomial-based  approximations, 
commonly  found  in  the  literature,  fail  for  pulsating  operation  con¬ 
ditions. 


2.  Single-particle  electrochemical  model 

The  considered  Li-Ion  cell,  presented  in  [19],  is  one  of  72  units 
taking  part  of  a  6  Ah  and  276  V  pack  designed  for  power  assistance 
in  hybrid  electric  vehicles.  Secondary  Li-Ion  cells  consist  of  three 
main  components:  negative  electrode,  separator  and  positive  elec¬ 
trode.  Electrodes  are  made  of  porous  materials  with  insertions  of 
lithium  and  inert  conductive  particles;  carbon  composites  are  used 
for  the  negative  and  metal  oxides  for  the  positive  electrode.  Separa¬ 
tor  is  also  a  porous  but  electronically  non-conductive  matrix  placed 
between  the  electrodes.  Pores  inside  the  three  elements  are  filled 
with  an  electrolyte  solution  based  on  a  lithium  salt,  setting  up  a 
ionic  path  all  along  the  cell. 

During  discharge,  lithium  ions  (Li+)  inside  the  active  material  of 
the  negative  electrode  diffuse  towards  the  boundary  between  solid 
and  solution  phases,  and  react  (deintercalation).  Once  in  the  elec¬ 
trolyte,  ions  travel  to  the  positive  electrode  where  a  corresponding 
reaction  takes  place  (intercalation),  but  now  ions  diffuse  from  the 
interphase  inside  the  solid  phase.  To  keep  an  electric  balance,  elec¬ 
trons  are  delivered  from  the  negative  electrode  and  received  in  the 
positive,  flowing  through  an  external  circuit  as  a  usable  current;  an 
inverse  and  complementary  process  occurs  when  charging. 

Electrochemical  models  of  Li-Ion  cells  are  based  on  the  porous 
electrode  theory  [36],  that  states  that  porous  electrodes  could  be 
treated  as  two  overlapped  phases.  Because  of  the  complexity  of 
pores  geometry,  distribution  of  variables  is  analyzed  just  in  the 
longitudinal  dimension  x  that  goes  along  the  three  elements  of 
the  cell.  Major  variables,  as  electrostatic  potential  in  both  phases, 
concentration  of  Li+  in  the  solid  matrix,  concentration  of  elec¬ 
trolyte  in  the  solution  and  ionic  flux  in  the  interphase  depend  on 
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both  time  and  position  x.  The  solid  phase  is  modeled  as  a  set  of 
regular  forms,  usually  spheres  of  the  same  radius.  Lithium  concen¬ 
tration  in  solid  phase  also  depends  on  r,  considered  as  a  second 
pseudo-dimension,  representing  the  position  along  the  radius  of  a 
representative  sphere  of  each  electrode. 

From  simulations  reported  in  [19,29,35]  of  the  cell  in  study,  it  is 
observed  that  ionic  flux  in  the  solid-solution  interphase,  as  well  as 
concentration  of  Li+  in  the  solid  phase,  become  nearly  independent 
of  position  x  along  electrodes  very  fast  after  the  beginning  of  dis¬ 
charge.  Furthermore,  if  these  two  variables  are  always  considered 
independent  of  position  x,  the  result  is  a  decoupled  model  where 
lithium  diffusion  in  the  solid  phase  (intercalation  and  deinterca¬ 
lation)  represents  the  dominant  dynamics  in  the  current-voltage 
relationship  of  the  cell.  In  this  case,  electrolyte  concentration  in 
the  solution  is  assumed  constant,  which  drives  to  a  simplification 
known  as  the  single-particle  electrochemical  model.  Such  model 
was  first  introduced  for  nickel-metal  hydride  (NiMH)  cells  [37], 
then  extended  for  Li-Ion  batteries  [26]  and  has  been  used  for  SOC 
estimation  in  works  as  [38-40].  Results  of  a  deduction  similar  to 
that  followed  in  [39]  are  summarized  below  and  a  more  formal 
approach  can  be  found  in  [20]. 

In  the  single  particle  model,  lithium-ion  diffusion  in  the  solid 
phase  is  studied  by  means  of  one  representative  particle  for  each 
electrode.  This  phenomenon  is  described  by  the  radial  diffusion 
equation 


—  -nl  — 

di~°^dr 
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losses  due  to  solution-phase  are  calculated  for  the  three  elements 
in  [39].  The  consecutive  addition  of  losses  results  in 
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where  Ssep  is  the  separator  thickness  and  /c,  with  the  corresponding 
subscript,  represents  the  effective  conductivity  of  the  electrolyte 
solution  in  each  element. 

In  [39],  overpotentials  pp  and  pn  are  found  from  the  non-linear 
Butler-Volmer  reaction  equation,  resulting  in  an  expression  equiv¬ 
alent  to 
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where  i0  is  the  exchange  current  density,  7 Z  the  ideal  gas  constant 
and  Tthe  absolute  temperature.  Alternatively,  using  the  linearized 
Butler-Volmer  equation 
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with  aa=ctc  =  0.5,  the  results  is 
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Taking  (9)  instead  of  (7),  voltage  at  terminals  is  finally 


(9) 


with  boundary  conditions 
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where  c  =  c(t,  r)  is  the  concentration  of  lithium  ions  inside  the  rep¬ 
resentative  particle  of  the  related  electrode,  D  is  the  diffusivity 
constant  of  Li+  in  the  solid  phase,  r  is  the  position  along  the  radius 
of  the  representative  particle  and  R  is  the  total  radius  of  the  parti¬ 
cle.  jn  =jn(t )  is  the  normal  ion  flux  in  the  solid-solution  boundary, 
averaged  along  the  associated  electrode  (see  [39])  such  that  it  is 
independent  of  position  x  and  proportional  to  the  applied  current 
/(f), 
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for  positive  and  negative  electrodes,  respectively.  8  is  the  electrode 
thickness,  A  is  the  sectional  area  of  the  cell,  T  is  the  Faraday’s  con¬ 
stant  and  a  is  the  electrode  specific  contact  area.  When  needed, 
subscripts  p,  n  and  sep  are  included  to  distinguish  between  positive, 
negative  electrodes  and  separator. 

Voltage  in  terminals  V  is  the  difference  of  solid  phase  potential 
0S  at  cell  edges: 


V  =  UP(0P)  —  Un(0n)  -  RqI, 


(10) 


with  the  total  resistance  Rq  defined  as 

R  =  *  (  sp  i  Ssep  [  Sn  nr  |  nr 
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In  summary,  the  single-particle  electrochemical  model  is  rep¬ 
resented  by  Eqs.  (l)-(3)  and  (10),  where  the  applied  current  I  is 
the  input,  the  voltage  V  at  terminals  is  the  output  signal  and  the 
state  variables  cn  and  cp  correspond  to  concentration  of  Li+  in  both 
electrodes.  Values  of  parameters,  extracted  from  [29],  are  listed 
in  Appendix  B.  For  pulsating  charge-discharge  conditions  around 
nominal  current,  assumption  of  uniform  electrolyte  concentration 
is  acceptable  in  general.  For  long  term  constant  charge  or  discharge 
regimes,  or  for  high  pulsating  currents,  it  would  be  necessary  to 
find  limiting  conditions  of  discharge  time  and/or  current  intensity 
to  assure  that  the  electrolyte  solution  along  the  cell  will  not  run  out 
of  electrolyte  in  any  place.  The  procedure  proposed  in  [25]  could 
be  followed  for  this  purpose. 

3.  Spatial  semidiscretization  approach 


)+«/•  ni) 


V  =  0s  lx=0  -  0s  \x=L  ~  Rfl,  (4) 

where  Rf  is  the  equivalent  film  resistance  at  solid-solution  inter¬ 
phase  of  electrodes  and  L  =  8P  +  8sep  +  <$n.  Averaging  the  cell  variables 
along  electrodes  and  separator,  Eq.  (4)  is  rewritten  as 

V  =  Up(6p)  —  Un(0n)  +  Pip  —  Pin  +  0e,p  —  0e,n  —  RfU  (5) 

with  overpotential  defined  as  p  =  0S  -  0e  -  U  and  0e  the  elec¬ 
trostatic  potential  in  the  solution.  Stoichiometry  0  =  c|r=p/cmax  is 
defined  as  the  normalized  concentration  of  Li+  in  the  solid-phase, 
also  known  as  the  electrode  state-of-charge  (do  not  confuse  with 
battery  state-of-charge  or  SOC).  Open  circuit  potential  functions 
1/(0)  of  both  electrodes,  taken  from  [19],  are  shown  in  Appendix  A 
as  well.  Usually,  the  process  of  electrolyte  solution  is  ignored;  in 
this  case  such  assumption  means  that  0e,p  -  0e,n  =  0.  Instead,  ohmic 


Spatial  semidiscretization,  also  known  as  the  method  of  lines, 
consists  on  approximating  the  state  variable  in  a  set  of  discrete 
points,  arbitrarily  located  in  the  spatial  domain,  by  means  of  a  sys¬ 
tem  of  ordinary  differential  and  algebraic  equations.  For  simplicity, 
the  state  and  position  variables  c  and  r  are  normalized  against  the 
maximum  value  cmax  and  the  radius  R  of  the  representative  particle, 
respectively.  Then  the  state  Eq.  (1 )  is  rewritten  as 
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and  boundary  conditions  (2)  take  the  form 


a{t)  = 


d© 

dp 


=  0  and  P{t)  = 


p= o 


d© 

dp 


p= i 


-Rjn 
D  ’ 


(13) 
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where  I<=DIR2,  ©(t,  p)  =  c(t,  r/ft)/cmax,  p  =  r/R  and  {©,  p}e[0,  1]. 
Electrode  state-of-charge,  which  is  needed  to  evaluate  the  objective 
function  (10),  is  then  0  =  ©(t,  1). 

Applying  the  spatial  semidiscretization  approach,  the  diffusion 
problem  (12)  and  (13)  is  approximated  for  each  electrode  by  the 
system  of  linear  ordinary  differential  equations 

^  =  K  (A©  +  B/3)  ,  (14) 

over  the  whole  set  of  points  pit  with  i  =  l,  2,  . . .,  N,  or  just  in 
a  subset.  Depending  on  the  discretization  method,  a  differential 
equation  does  not  exist  at  specific  points,  usually  the  bound¬ 
aries.  This  happens  for  the  finite  difference  and  the  differential 
quadrature  methods  applied  in  a  classical  fashion.  Electrode  state- 
of-charge  0  =  0(t,  1 )  =  ©jv,  which  lies  in  the  outermost  boundary,  it 
is  calculated  as  a  linear  combination  of  the  elements  of  G  and  the 
boundary  condition  /?,  such  that 

0  =  C0  +  D/3.  (15) 


The  value  of  state  variable  at  p  =  0  is  not  of  interest  for  the  single¬ 
particle  model.  However,  ©(t,  0)  =  ©!  is  needed  for  evaluations  in 
time  domain  and  could  be  calculated  similarly  to  ©^. 

Although  0jvG0  for  the  finite  element  method,  representation 
(1 5)  is  also  valid  for  such  discretization  method.  Matrices  A  e  Rqxq} 
B  e  Mqxl ,  C  e  M1  xq  and  DeRarea  result  of  discretization  and  their 
elements  are  different  in  general,  depending  on  the  related  elec¬ 
trode,  the  discretization  method  and  the  number  of  points  as  well 
as  their  distribution.  Matrices  B  and  D  are  the  second  column  of 
BeRqx2  and  DeMlx2,  respectively.  The  first  column  of  B  and  D 
corresponds  to  the  boundary  condition  a ,  which  is  out  of  the  for¬ 
mulation  due  to  its  null  influence.  Index  q,  defined  later  for  each 
discretization  method,  is  the  length  of  vector  G  and  also  represents 
the  order  of  the  system  of  differential  Eq.  (14).1  Thus,  qp  is  the  order 
of  the  subsystem  which  approximates  the  diffusion  equation  asso¬ 
ciated  to  the  positive  electrode  and  qn  is  similarly  defined  for  the 
negative  electrode. 

According  to  the  spatial  semidiscretization  approach,  the  com¬ 
plete  single-particle  electrochemical  model  is  approximated  by  the 
system  of  ordinary  differential  and  algebraic  equations 


©p] 

KpAp  0 

QP 

—jipBp 

©n 

0  KnAn 

Qn 

\ 

P'nBn 

0p  =  CpGp  +  DpI , 

6>„  =  Cn0n  +  Dn/, 

V  =  Up(0p)  -  Un(0n)  -  RqI, 


(16) 

(17) 

(18) 
(19) 


where  /x  =  R/iDa^SA).  In  (16),  dynamics  of  positive  and  negative 
electrodes  are  put  together  and  then  the  order  of  the  complete 
approximated  model  is  Q  =  qp  +  qn- 


4.  Discretization  methods 


problem  (12)  and  (13):  finite  difference,  finite  element  and  differ¬ 
ential  quadrature.  In  addition,  two  criteria  for  order  selection  are 
explored  and  compared.  In  these  subsections,  a  brief  overview  of 
the  considered  discretization  methods  is  presented.  All  of  them  are 
applied  in  a  classical  fashion. 


4A.  Finite  difference  method  (FD) 


The  three-point  centered  finite  difference  is  maybe  the  most 
popular  discretization  method.  In  this  case,  for  one  dimension  prob¬ 
lems,  derivatives  of  any  order  with  respect  to  x  of  a  function  f[x) 
are  approximated  throughout  a  grid  by  means  of  point  operators, 
which  depend  on  f[x)  evaluated  in  the  point  of  interest  x2  and  the 
neighboring  points  xz_!  and  x2+1 .  Operators  are  derived  by  taking 
linear  combinations  of  Taylor  series  of/(x2_i  ),/(x2)  and/(x2+1 )  around 
x2,as  well  as  3/(xi  )/3xand  df[xN)ldx  when  necessary.  See  [41,42]  for 
more  details  and  further  FD  schemes. 

With  the  three-point  centered  FD  method  it  is  not  possible  to 
formulate  differential  equations  at  the  boundary  points  p\  =  0  and 

Pm  =  1,  then,  state  vector  is  0  =  [©2  ©3  -  -©jv-i]  with  ©2  =  ©|Pi. 
and  q  =  N  -  2.  It  is  convenient  to  take  the  expanded  form 


3©  /a2©  23©\ 

3 1  ~  \dp2  +  p  dp) 


(20) 


of  Eq.  ( 1 2 ),  whose  discretization  with  the  introduction  of  boundary 
conditions  results  in 


d©2- 

df 


=  I<  < 


2  +  1 

E 

I,  z 

1  ®k  +  ^ct 

for 

i  =  2, 

K=l 

2+1 

E  + 

t  +  —  ai,l< 
Hi 

)©k 

for 

z  =  3,..., 

li=2-l 

t  + 

k=i- 1 

2 

^  +  ~77ai,k 
Pi 

)  ®k  +  CnP 

for 

i  =  N-  1, 

(21) 

where  N  is  the  total  number  of  points.  State  variable  ©(t,  p)  is 
approximated  at  boundaries  p  =  0  and  p  =  1  through  the  expressions 


©i  =  d\oc  +  d2©  2  +  d3©3, 


(22) 


©jv  =  djv-2©iv-2  +  djv_i  ©jv-i  +  djv/3,  (23) 

that  can  be  obtained  from  a  combination  of  the  first  and  second,  and 
second  and  third  lines  of  (21 ),  respectively.  Coefficients  aik  and  bik 
are  related  to  the  first  and  second  spatial  derivative  and,  in  addition 
to  Qv  and  d0,  are  explicitly  calculated  in  function  of  discretization 
points  as  shown  in  Table  1,  whereas  C\  is  irrelevant  here  and  in  all 
cases  as  a  =  0.  Once  coefficients  are  known,  (21)  and  (22)  could  be 
easily  brought  to  matrix  form  (14)  to  (15). 


Some  approximations  based  on  spatial  semidiscretization  for 
diffusion  problem  (1)  and  (2)  have  been  reported.  For  example,  a 
fifth  order  finite  element  approximation  is  proposed  in  [35]  and 
the  finite  difference  method  is  applied  in  [39].  Nevertheless,  in  nei¬ 
ther  of  both  works  it  is  specified  why  to  choose  one  or  the  other 
discretization  method  or  how  to  determine  the  order  of  approxi¬ 
mation.  Thus,  the  aim  of  this  work  is  to  identify  which  one  among 
three  discretization  methods,  combined  with  four  grids,  offers  the 
best  trade  off  between  low  order  and  accuracy  for  the  diffusion 


1  Unless  stated  otherwise,  in  this  paper  order  always  refers  to  the  number  of 

differential  equations  in  the  resulting  system. 


4.2.  Finite  element  method  (FE) 

This  is  also  a  well  known  discretization  method  and  is  widely 
used  to  simulate  multidimensional  problems  with  irregular  geome¬ 
tries  because  of  its  easy  numerical  assembling.  For  one  dimensional 
problems  and  taking  the  simplest  formulation,  the  aim  of  the  FE 
method  is  to  approximate  the  profile  of f[x)  by  linear  segments  like 

m  «  +  -^3-/(xi+ 1 ),  (24) 

*2+1  *2  *2+1  *2 

where  the  factors  of/(x2)  and/(x2+i )  are  called  aspect  functions.  With 
this  method  it  is  possible  to  define  ordinary  differential  equations 
in  every  point  when  Neumann  type  boundary  conditions  as  (13) 
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Table  1  4.3.  Differential  quadrature  method  (DQJ 

Coefficients  for  the  FD  method. 


i  Coefficients 


1 

2 

3 


N-2 


N-  1 


N 


d,  =  _ ^3 -Pi)2 _ 

(P3  ~P2  )(P3  +P2  -2Pl ) 

^2’2  (P3-P2)(P3+P2-2Pl ) 

b2,3  =  -  b2,2 

U.  .  _  _ 2 _ 

1,1-1  (Pf-Pi-1  XPj+1  -Pi-1  ) 


1,1  (Pt— Pj-lXPi+l-Pj) 

^M+1  =  (Pi+l-Pi)(Pi+i-Pi-i) 


A  _  ~ (P2  ~ Pi  f 

J  (P3-P1XP3+P2) 

=  _ ~ 2(P2— Pi ) _ 

(P3-P2)(P3+P2-2Pi) 
02,3  =  -  <32,2 

_  -(Pi+l-Pj) 
Ul’1-1  “  (Pf-Pi-1  XPi+i -Pi-i ) 

_  Pj+i  — Zpj+Pj-! 
UM-1  -  (Pi-p1-_1  XPi+1  -Pi) 

_  Pj~  Pj- 1 

01,1+1  -  (Pi+1-Pi)(Pi+i-Pi_i) 


hiV  1,N  2  (PiV-1  — PN— 2  )(2pj\f— PN-1  —PN—2 ) 


n  _  _ —  2(pjy—  PN—  l) _ 

N—l,N—2  (pN_1  — PjV— 2  )(2pj\f— PiV-1  -PN-2) 

friV-l,N-l  =  -  bN-l,iV-2  <3N-1,N-1  =  ~  ClN-l,N-2 


c  _  2+pN_T  -Pjy_2 
N  2PN-PN-i-PN-2 

d  _  _ -(PW-PAf- 1)2 _ 

™  Z  (Pn-1 -PN-2)(2Pn-PN-1-PN-2) 

,  =  _ (PN-PN-2)2 _ 

™  1  (PN-1-PN-2)(2Pn-Pn_i-Pjv_2) 

j  _  (PN-PN-l)(PN-PN-2) 

N  (2Pn-PN-1-PN-2) 


are  given.  In  this  case  0  =  [©1  ©2  •  •  -©n]  T  and  q  =  N.  The  diffusion 
problem  (12)  and  (13)  is  approximated  by 


r  i+1 

Ebi-kQfc+clg 

for 

/  =  1, 

k=i 

2+1 

k=i-\ 

for 

i  =  2, . . . ,  N  —  1, 

(25) 

^  +  CN^ 

<  k=i-\ 

for 

i  =  N-  1. 

In  contrast  to  both  methods  described  before,  differential 
quadrature  is  less  popular  despite  its  generality,  easy  implemen¬ 
tation  and  good  approximations.  It  is  based  in  the  assumption 
that  some  smooth  function  /(x)  can  be  approximated  by  a  poly¬ 
nomial  like  X^=i^xk_1  or  a  trigonometric  series  of  form  A0  + 

Ylk=i^k  cos  k*  +  Bk  sinfcx).  The  first  case  is  called  polynomial  dif¬ 
ferential  quadrature  (PDQ)  and  the  second  case  corresponds  to  the 
Fourier-based  differential  quadrature  (FDQ).  As  a  consequence  of 
this  formulation,  any  differential  operator  dnf[x)ldxn  of  order  n  is 
approximated  at  xz  by  linear  combinations  of  the  functional  values 
of /(x)  at  every  discretization  point  [44],  not  just  at  adjacent  points 
as  for  FD  and  FE. 

When  applied  in  a  classical  way,  no  differential  equations  can 
be  stated  for  the  boundary  points  with  the  DQ method  (see  [45]  for 
imposition  of  boundary  conditions  at  interior  points),  and  (9  and  q 
are  the  same  that  for  FD.  In  general,  the  expanded  form  (20)  of  (12) 
is  discretized  as 


d©t-  _ 
dr  - 


~N-1 

^  ^  (^i,k  +  ~^i,k^j  ®k  +  CiP 


for  i  =  2,  ...,N  — 1. 

(27) 


It  should  be  underlined  that  there  is  a  coefficient  q  for  each  point 
because  ft  has  a  direct  influence  in  every  location. 

Due  to  imposition  of  boundary  conditions,  coefficients  akk,  bik 
and  Q  are  indirectly  calculated  in  function  of  the  discretization 
points  Pi  through  the  more  elemental  coefficients  aik  and  bik, 
related  to  the  first  and  second  derivative  as  well,  such  that 

bi,k  =bi,k+P  [al,k(bi,NaJV,l  -  &i,lflN,N)  +  «N,k(^i,lal,N  -  ^i,Nal,l)]  > 

(28) 


Coefficients  bi<k  and  Qv,  shown  in  Table  2,  are  found  by  solving  the 
error  integral  function 


Ri  =  - 


Wi 


P 2  3  P 


9© 

dt 


(26) 


for  each  grid  point,  including  the  boundaries.  The  integral  is  evalu¬ 
ated  in  two  parts,  one  for  each  element  associated  to  pit  and  taking 
the  appropriate  aspect  function  as  weight  Wj-  (Galerkin’s  method). 
For  the  temporal  derivative  of  ©,  the  lumped  formulation  has  been 
considered,  which  consists  on  assuming  the  term  9©/9t  as  a  con¬ 
stant  for  the  point  p\  of  interest  and  half  of  each  of  the  two  adjacent 
elements  see  [43].  The  consistent  formulation  has  been  discarded 
because  of  its  probed  tendency  to  produce  oscillations. 


Table  2 

Coefficients  for  the  FE  method. 


fii,k  =  ai,k  +  P  [al,k(ai,NaN,l  -  Oi.lflN.fv)  +  dN,kiai,\a\,N  ~  1i,Na^P  )]  - 

(29) 


Ci=P 


bi,Nal,1  -<)i,lal,N+  —(di,Na  1,1  -  fy, iQi.n) 

Pi 


(30) 


where  P=(0i,i0n,n-0i,n0n,i)_1-  At  boundaries,  state  variable  ©(t, 
p)  is  approximated  by  the  equations 


©!=P 


N-l 

,nP  +  ^^(01  ,N0JV,k  -  0JV,N01,k)©k 

k= 2 


(31) 


®N  =  P 


N- 1 

-On, \P  +  y^(0N,10i,k  -  01,1 0N,k)®/< 
k= 2 


(32) 


N-l 

N 


Coefficients 


&i,i  =  ■ 


hi- 1  = 


hi  = 


|  \{P\+P2)3-P\^  (Pl-P2)2 
-8(Pi3_1-P^) 

[-(Pi_i  +Pf  )3  +(Pj -Pj-i  )3  J  (Pi_ i  +Pi )' 


hjv.N-i  = 


-(Pi_i+P/r+(p,—Pi+i r  I  (Pf_i  +Pj)  . 
[-(Pi-i  +Pi)3+(p,-Pi+i  )3  j  (Pi+Pi+i  )2 


hi,2  =  -  hi, 


AizL  +  _++■, 

(Pi_i  —Pi)2  (Pj-Pi+iX 


[44(^v- 


i-Pn)  (pn-i-pn) 


h/v,N  —  hjvjv-i 


Qv  =  -f 


5  ^N-l  +PN  )3 


In  the  case  of  PDQ,  coefficients  az  k  are 


0;,k  = 


fj  (Pi-Pj) 

J=1  ,j  +  i 

N 

(. Pi-Pk )  fj  ( Pk-Pj ) 


-  E  a‘-fe 

k=l,k±i 


for  i^k,  {i,k}  =  i,...,N 

for  i  =  k,  {i,k}  =  1, ...  ,N. 

(33) 


Likewise,  for  FDQ 
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N 

]  J  sin ((a-Pj)/2) 
J=l 


di,k  =  < 


sin((p/  -  pk)/2)  JJ  sin((pk  -  pf)/ 2) 
j— i  J  ^  ^ 

N 

-  iz  a,-k 


V  k=1,k^i 


for 


for  i  =  k. 


{i,k}  =  1,  ...,N 


{/,/<}  =  1,  ...,N. 

(34) 


In  both  cases,  if  coefficients  a2i/<  are  the  elements  of  a  matrix  W(1), 
then  bi  k  is  an  element  of  matrix  W(2)  = 


Fig.  1.  Energy  spectrum  of  a  rectangular  pulse. 


4.4.  Distributions  of  discretization  points 


5.1.  Truncation  by  bandwidth  selection 


Accuracy  of  discretization  methods  depends  on  a  great  extent 
on  the  location  of  discretization  points.  For  this  reason,  four  grids 
are  tested  along  with  the  discretization  methods. 

Uniform  distribution  (UD):  Discretization  points  are  evenly 
placed  throughout  the  domain  of  interest. 

Chebyshev-Gauss-Lobatto  distribution  (CGL):  points  become 
closer  towards  the  domain  boundaries  according  to 

Pi  =  1  (l  -  COS  S-f for  i  =  l,...,N.  (35) 

Roots  of  shifted  Legendre  polynomials  (RLP):  Legendre  associated 
polynomials  are  deducted  from  the  recursive  formula 

p„+m  =  (2n+1)^( p)-”Pn-mt  (36) 

with  P0  =  1  and  ?!  =  p.  Roots  of  such  polynomials  lie  on  p  e  (-1 , 1 ). 
After  shifting  the  roots  to  the  domain  pe(0,  1)  by  means  of  the 
transformation  p  =  2p  -  1 ,  and  including  the  auxiliary  points  p\  =  0 
and  pN  =  1,  the  elements  of  the  resulting  set  could  be  used  as  dis¬ 
cretization  points. 

Proportional  distribution  (PD):  As  stated  in  [43],  results  from  FE 
method  can  be  improved  if  the  grid  is  tightened  where  changes  in 
the  variable  of  interest  are  more  pronounced.  For  this  grid,  points 
become  closer  towards  the  boundary  point  p  =  1  in  a  square  pro¬ 
portion  defined  by 

«  -  \fff-  <37) 

Such  heuristic  distribution  is  proposed  because  it  is  known  that 
the  state  variable  tends  to  have  a  parabolic  profile  in  the  spatial 
variable. 

5.  Order  selection 


This  criterion  consists  on  choosing  the  order  of  the  system  such 
that  its  spectrum  covers  all  or,  at  least,  the  major  part  of  the  band¬ 
width  associated  to  a  test  signal.  Supposing  (5  is  a  rectangular  pulse 


«')-{; 


for 

for 


te[ 0,  d) 

t  e  [ d ,  oo), 


(38) 


of  unitary  amplitude  and  duration  d  >  0,  it  can  be  easily  shown  that 
its  energy  spectrum  is  given  by 


E[o>)  = 


=  d2sinc2 


(39) 


where  is  the  Fourier  transformation  of  (38)  and 

sinc(x)  =  sin  (ttxJ/ttx.  From  Fig.  1  it  is  observed  that  the  first 
zero  of  (39)  occurs  at  (jo  =  ±2jrld.  Furthermore,  taking  the  integral 
of  (39)  over  the  interval  co  e[-2nld,  2jr/d],  it  is  concluded  that 
~90%  of  energy  is  contained  in  the  bandwidth  oo  e  [0,  2jr/d]. 

The  characteristic  values  km  of  the  diffusion  problem  (12)  and 
(13)  could  be  found  solving  the  associated  regular  Sturm-Liouville 
problem  for  the  homogeneous  case  [46].  All  are  real  and  nonnega¬ 
tive  [35]  such  that 


f  0  for  m  =  1 

\  ym  >  0|ym  =  tanym  for  m  >  1 , 


(40) 


and  the  transient  response  of  the  solution  @(t,  p)  of  (12)  and  (13) 
to  any  input  will  be  determined  by  modes  of  the  form  exp  (pmt), 
with  pm  =  for  m  =  1 , 2, . . .  Let  pq  be  the  constant  of  the  faster 

mode.  Then,  with  this  truncation  criterion,  it  should  be  satisfied 
that  | pq  |  >  2nld  to  ensure  that  the  reduced  order  model  spectrum 
covers  the  bandwidth  where  most  of  energy  of  the  input  signal  is 
concentrated,  which  implies  that 


(41) 


The  two  order  selection  or  truncation  criteria  suggested  but  not 
discussed  in  [30]  are  explored  in  this  section.  The  first  is  based  in 
the  model  bandwidth  selection  whereas,  in  the  second,  a  residue 
analysis  is  carried  out.  In  both  cases,  the  supposition  of  a  specific 
profile  for  the  boundary  condition  /3(t),  the  input  of  each  electrode 
subsystem,  plays  a  key  role.  In  particular,  it  is  considered  that  the 
current  profile  used  in  [35],  the  FreedomCAR  Flybrid  Pulse  Power 
Characterization  (HPPC)  [34],  is  applied  to  the  Li-Ion  cell.  One  cycle 
of  that  test  signal  consists  on  three  constant  stages:  a  30  A  and  18  s 
discharge  pulse,  an  open-circuit  (0  A)  relaxation  period  of  32  s  and, 
finally,  a  22.5  A  and  10s  charge  pulse.  Nevertheless,  it  is  helpful  first 
to  study  the  case  where  is  a  simple  rectangular  pulse  as  results  are 
straightforward  extended  for  every  signal  composed  of  rectangular 
pulses  of  any  amplitude  and  duration. 


When  p  consists  of  more  than  one  pulse,  it  is  enough  just  to  study 
that  of  minimum  duration  because  its  90%-energy  bandwidth  actu¬ 
ally  covers  more  than  that  for  longer  duration  pulses. 

5.2.  Truncation  by  residue  analysis 

The  output  function  (10)  or  (19)  depends  on  the  concentration 
of  Li+  evaluated  just  at  the  surface  of  the  representative  particle 
of  each  electrode,  namely,  ?  =  ©(t,  1).  If  variables  ft  and  ©(t,  1) 
are  defined  as  input  and  output  signals,  respectively,  the  transfer 
function  (TF)  which  defines  the  response  of  ©(t,  1)  to  ft  is 

G(s,l)=gfrl)=  -tanlKJfc)  ^  (42) 

P(s)  tanh(VKs)  -  VKs 
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where  fi(s)  and  0(s,  1)  are  the  Laplace  transformations  of  fi(t)  and 
@(t,  1)  (see  [30,35]). 

Transfer  function  G(s,  1)  is  also  understood  as  the  impulse 
response  of  0(s,  1 )  and  its  modal  expansion  results  in  the  infinite 
series 


G(s,  \  )  =  K 


2Z 


1 


S-Pm 


(43) 


whose  poles  pm  are  defined  above.  This  expression  shows  that 
modes  related  to  nonzero  poles  are  equally  important  and  no  trun¬ 
cation  criteria  could  be  established  from  the  impulse  response. 
Applying  the  input  pulse  (38),  with  Laplace  transformation 

jS(s)  =1(1—  e“ds)  ,  (44) 

the  expansion  of  the  convolution  between  C(s,  1)  and  fi(s)  is 


G(s,  1  )jS(s)  =  K 


m= 2 


1 

S-Pm 


(45) 


In  this  case,  each  mode  is  weighted  by  a  residue  resm  =  2/C/pm,  which 
decays  as  its  associated  mode  becomes  faster.  This  fact  allows 
to  define  a  truncation  criteria.  Clearly,  the  most  important  mode 
related  to  a  nonzero  pole  is  the  second  (m  =  2).  Then,  the  order  of 
the  approximation  in  (14)  could  be  selected  such  that 


resq 

P2 

res2 

Pq 

(46) 


6.  Other  approximations 

Whereas  spatial  semidiscretization  is  a  very  general  approach, 
most  of  reduced  order  approximations  reported  in  the  literature 
for  the  solid  phase  diffusion  problem  (1)  and  (2)  are  valid  just 
under  constant  charge  or  discharge  input  currents.  For  instance, 
an  expression  based  on  the  Dunhamel’s  superposition  integral, 
valid  either  for  short  or  long  periods  of  operation,  is  proposed 
in  [16,47]  to  obtain  more  efficiently  numerical  solutions  of  the 
whole  electrochemical  model.  Simple  approximations  adjusted  to 
the  analytic  solution  are  presented  in  [35,48],  also  for  short  time 
and  steady-state  predictions.  In  [17,25]  an  equation  representing 
the  steady-state  solution  is  used  to  analyze  the  solid  phase  diffusion 
limitations  and  determine  the  capacity  of  a  cell. 

For  time  varying  input  rates,  two  dynamic  approximations 
based  on  low  order  polynomials  are  presented  and  explicitly  solved 
for  some  input  profiles  in  [27].  For  the  first,  it  is  assumed  that 
lithium  concentration  c  can  be  described  by  a  second  order  poly¬ 
nomial  in  r  of  form  c(r,  t)  =  a(t)  +  b(t)(r/ft)2.  Solving  for  weights  a(t) 
and  b(t)  results 

dc  -3 . 

df  =  ~R  h  R  (50) 

C|r=R  =  ~C  ~  5 Djn’ 

which  is  a  first  order  approximation  and  just  keeps  the  dynamics 
associated  to  the  average  lithium  concentration  c  inside  a  spheric 
particle. 

Taking  the  fourth  order  polynomial  c(r, 
t)  =  a(t)  +  b(t)(r/ft)2 +  d(t)(r/ft)4  and  now  solving  for  a(t),  b(t) 
and  d(t)  results  the  second  order  approximation 


where  e  <  1  is  an  arbitrary  parameter  that  could  be  interpreted  as 
an  error  measurement. 

Moreover,  if  an  oscillating  input  fd{t)  =  sin  ( (Dot )  is  chosen,  the 
expanded  response  results  as 


C(s,  1  )/3(s)  =  K  [  — +2V- 

\  (DoS  2L^( 


(Dq 


(DoS 

where  coefficients 

(Dq 


®0  +  Pm  S~Pm  S2  +  ®o 


1  b0S  +  C0 


171=2 


b0  =  —  -  2  V 

(D0  ^ 


m= 2 


col  +  P™ 


and  c0  =  -2^~^ 

m=2 


Pm 


COq  +  Pm 


(47) 


(48) 


define  the  steady-state  response.  In  this  case,  residues  resm  = 
2K(d0/{(d2  +  ph)  are  also  functions  of  frequency  <d0  and  take  their 
maximum  value  resm  =  2K/\pm\  at  oj0  =\pm  I-  Then,  criterion  (46) 
can  be  formulated  taking  the  maximum  residues  resm: 


res  q 

P2 

res2 

Pq 

(49) 


In  this  case,  the  response  is  negligible  whether  cu0>  I  Pq  |.  When  the 
pulse  (38)  is  considered  as  an  infinite  sum  of  oscillating  functions, 
interpretation  of  this  last  result  could  be  that  the  response  due  to 
components  of  frequency  higher  than  \pq\  could  be  neglected. 

Based  on  a  modal  grouping  of  an  expansion  similar  to  (45), 
obtained  for  a  unit  step  input  signal,  a  low  order  and  large  band¬ 
width  transfer  function  approximated  to  (42)  was  synthesized  in 
[30].  In  contrast,  expansions  (43),  (45)  and  (47)  are  used  in  this 
paper  to  show  that  simple  residue  analysis  could  brought  a  prac¬ 
tical  order  selection  criterion  for  the  discretized  state  equations. 
It  should  be  shown  later  that  reasonably  results  are  straightfor¬ 
ward  obtained  without  further  alterations  of  the  semidiscretization 
approximations,  as  could  be  suggested  by  the  residue  grouping 
method. 


dc  _  -3  . 

d t  =  ~RJn 

dq  _  -30 D  _  45  . 

dt  -  R2  Q  ~  2R2Jn 
8  R  .  R 

C\r=R  =  C-  ™  <7-  or! 


(51) 


where  q  is  the  volume-averaged  concentration  flux,  which  defines 
the  average  change  of  concentration  c(t,  r)  with  respect  to  the  posi¬ 
tion  r.  q  is  calculated  by  integrating  of  3 c/3r  along  r  on  the  domain 
re[ 0,  R]  (see  [27]  for  more  details).  This  approach  is  applied  to 
the  single-particle  model  for  estimation  of  the  state-of-charge  in 
[38,40],  state-of-health  in  [28]  and  for  validation  of  reduced  order 
models  through  numerical  simulations  in  [20,31,33].  Because  of 
their  recurrent  appearance  in  the  literature,  (50)  and  (51)  are  also 
tested  and  compared  with  the  semidiscretization  approximations 
described  above. 

Some  other  approximations  have  been  proposed  in  the  fre¬ 
quency  domain.  In  [29]  the  exact  transfer  function  from  Jn(s)  to  C(s, 
R),  the  ratio  of  the  Laplace  transformations  of jn(t)  and  c(t,  R)  =  c\  r=R, 
is  approximated  through  a  lumped  transfer  function  whose  order  is 
arbitrarily  proposed  with  poles  optimally  found  for  a  certain  band¬ 
width.  A  similar  approach  is  developed  in  [30]  although,  in  this  case, 
poles  are  found  by  grouping  in  frequency  belts  the  residues  associ¬ 
ated  to  a  unit  step  input  and  taking  the  weighted  average  of  their 
amplitude. 


7.  Results  and  discussion 


7.1.  Order  selection 


Truncation  criteria  brought  contrasting  results.  The  first  cri¬ 
terion  is  too  conservative  and  drives  to  relative  high  order 
approximations;  in  this  case  order  means  the  number  of  ordinary 
differential  equations  resulting  in  the  approximations.  For  the  neg¬ 
ative  electrode  results  that,  according  to  (41),  X2  q  >  3 141. 6  and  the 


10274 


A.  Romero -Becerril,  L.  Alvarez-Icaza  /  Journal  of  Power  Sources  196  (201 1 )  10267-10279 


first  value  which  accomplish  that  condition  is  X2  19  =  3375.9.  Simi¬ 
lar  results  are  obtained  for  the  positive  electrode.  In  this  case  X2q  > 
1696.5  and  then  X2  14  =  1796.7.  In  summary,  the  negative  elec¬ 
trode  subsystem  results  of  order  qn  =  19  whereas,  for  the  positive 
electrode,  the  subsystem  is  of  order  qp  =  14.  The  whole  approxi¬ 
mated  single-particle  electrochemical  model  is  of  order  Q  =  33. 

The  plot  of  normalized  residues  |resp?m/resp?2l  =  ^2 Am  °f 
expansion  (45)  against  the  index  m  is  shown  in  Fig.  5.  Defining 
arbitrarily  £  =  0.05  in  (46)  means  that  amplitude  of  the  residue  asso¬ 
ciated  to  the  faster  characteristic  value  will  be  5%  the  amplitude  of 
that  for  the  slower  mode.  Taking  such  value  for  £  the  result  is  that 
subsystems  of  both  electrodes  would  be  of  order  qn  =  qp  =  7  because 
X^/X2  =  X^/Xj  <  0.05  is  held  as  required  in  condition  (46).  Then, 
the  whole  approximation  is  of  order  Q  =  14,  which  is  less  than  half 
of  the  order  calculated  with  the  bandwidth  criterion. 

Approximations  are  first  evaluated  in  the  time  domain  to  have 
a  preliminary  surmise  about  which  of  them  has  the  best  trade  off 
between  low  order  and  accuracy.  In  this  case,  the  test  signal  con¬ 
sists  in  a  rectangular  pulse  followed  by  a  relaxation  time.  Then, 
the  frequency  response  of  approximations  is  compared  against  that 
of  the  analytic  transfer  function  (TF)  related  to  the  state  equation 
in  order  to  confirm  the  results  obtained  in  the  time  domain  tests. 
Finally,  representative  cases  of  all  considered  approximations  are 
applied  to  the  single-particle  model  and  their  response  to  the  FIPPC 


test  signal  is  compared  against  experimental  data  extracted  from 
[35]. 


7.2.  Time  domain  tests 

Time  was  normalized  such  that  r  =  I<t,  where  I<=D/R2,  and 
experiments  were  carried  out  in  the  interval  re[0,  0.45].  The 
pulse 


fl  for  re [0,0.15) 

P[  J  \0  for  re  [0.15,  0.45] 


(52) 


was  selected  as  test  signal  because  it  allows  the  system  to  reach 
the  steady-state  when  a  constant  input  is  forced,  /3(  r )  =  1 ,  as  well  as 
during  the  relaxation  stage,  /J(r)  =  0.  Programming  was  carried  out 
with  MathWorks®  MATLAB®.  The  base  solution  was  numerically 
obtained  through  the  pdepe  routine,  which  is  an  initial-boundary 
value  problem  solver  for  parabolic-elliptic  partial  differential  equa¬ 
tions  (PDEs)  in  one  dimension,  based  on  the  algorithm  proposed 
in  [49].  The  algorithm  employed  by  the  solver  performs  a  spatial 
semidiscretization  of  parabolic  and  hyperbolic  PDEs  by  a  piecewise 
nonlinear  Galerkin/Petrov-Galerkin  method  whose  truncation 


a  Total  error  for  finite  difference  (FD)  b  Total  error  for  finite  element  (FE) 


Approximation  order  q 


c  Total  error  for  polynomial  d  Total  error  for  Fourier-based 

differential  quadrature  (PDQ)  differential  quadrature  (FDQ) 


Fig.  2.  Results  of  tests  in  the  time  domain. 
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a  Time  profiles  for  finite  difference  (FD) 


b  Time  profiles  for  finite  element  (FE) 


C  Time  profiles  for  differential  quadrature  (PDQ) 


d  Time  profiles  for  Fourier-based  differential  quadrature  (FDQ) 


Fig.  3.  Time  profiles  for  approximations  of  order  q  =  7  using  proportional  distribution  (PD). 


error  is  of  second  order ?  A  vector  of  200  elements  uniformly  dis¬ 
tributed  on  the  space  domain  was  chosen.  Time  integration,  for 
the  base  solution  and  all  the  approximations,  was  solved  with  the 
odel5s  routine,  a  variable  order  and  multistep  implicit  algorithm 
for  stiff  problems,  which  is  based  on  numerical  differentiation  (see 
[50]  for  more  details).  Optionally,  the  solver  uses  the  backward  dif¬ 
ferentiation  formulas  known  as  Gear’s  method.  A  time  vector  of 
1000  elements  evenly  separated  was  specified  for  the  base  solu¬ 
tion  and  the  approximations.  The  solver  includes  a  scalar  relative 
error  tolerance  and  a  vector  of  absolute  error  tolerances  for  which 
the  default  values  1  x  10-3  and  1  x  10-6  were  kept,  respectively. 

Each  discretization  method  was  tested  along  with  the  four  point 
distributions  described  in  Section  4.4.  The  total  error  was  individ¬ 
ually  calculated  for  each  combination  using  the  function 


Et  = 


£2{€>r(0  -  @base(0 1 
^2{®base(0} 


x  100%, 


(53) 


where  vector  Ot(  t)  is  formed  by  the  state  vector  6>(  t)  and  the  values 
of  ©i(t)  and  ©N(t)  related  to  the  currently  tested  approximation, 
@base  results  from  interpolating  the  base  solution  on  the  corre¬ 
sponding  grid  and  £2!°}  implies  the  Lebesgue  2-norm.  Interpolation 
was  carried  out  with  the  pdeval  routine,  which  approximates  the 
variable  0(t,  p)  and  its  derivative  9©(t,  p)/dp  at  any  specified  point 
using  the  pdepe  routine. 


2  In  this  case  order  is  related  to  the  truncation  error  of  the  numerical  method. 


Error  of  approximations  evaluated  in  the  time  domain  are 
shown  in  Fig.  2.  Excluding  the  FE  method,  approximations  show 
their  best  performance  when  RLP  were  used,  but  similar  results  are 
observed  for  the  CGL  distribution,  which  is  easier  to  calculate.  In 
particular,  accuracy  of  PDQ  method  stands  out,  even  for  low  order 
approximations  as  shown  in  Fig.  2c.  This  result  is  expected  as  ordi¬ 
nary  differential  equations  resulting  from  applying  the  differential 
quadrature  method  are  strongly  coupled.  In  other  words,  informa¬ 
tion  from  every  discretization  point  is  used  to  infer  the  value  of  the 
state  variable  at  any  point  of  interest.  For  approximations  (50)  and 
(51 )  based  in  2nd  and  4th  order  polynomials,  because  of  no  grid  is 
related  to  these  approaches  and  concentration  of  Li+  is  of  interest 
just  at  the  surface  of  representative  particles,  error  was  calculated 
just  at  p  =  1  resulting  of  9.12%  and  1.57%,  respectively.  For  the  first 
case,  error  is  too  high  whereas  the  second  approximation  is  better 
than  any  of  the  same  order. 

In  general,  for  all  approximations,  error  ET  decreases  whereas 
order  q  increases  as  expected  for  consistent  methods.  An  excep¬ 
tion  occurs  for  PDQ  and  FDQ  when  PD  is  used;  in  this  case, 
those  methods  have  their  worst  error.  PD  produces  large  errors 
in  p  =  0  an  neighboring  points  as  shown  in  Fig.  3c  and  d.  This 
phenomenon  could  be  explained  because  points  become  more  sep¬ 
arated  and  with  similar  ordinate  towards  p  =  0.  Approximation  of 
spatial  derivatives  of  boundary  conditions  and  state  equation  lose 
accuracy,  inducing  growing  errors  for  high  order  approximations 
as  shown  in  Fig.  2c  and  d.  Whether  prior  knowledge  about  the  spa¬ 
tial  profile  of  state  equation  exist,  it  should  be  a  good  choice  to 
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a  Magnitude  frequency  response  of 
worst  approximations 


b  Phase  frequency  response  of 
worst  approximations 


c  Magnitude  frequency  response  of 
best  approximations 


d  Phase  frequency  response  of 
best  approximations 


Fig.  4.  Results  of  tests  in  the  frequency  domain. 


tighten  up  the  grid  where  more  pronounced  changes  are  found  as 
suggested  in  [43].  Thus,  PD  results  to  be  favorable  for  FE  as  shown 
in  Figs.  2b  and  3b  but  does  not  benefit  any  other  method.  On  the 
other  hand,  the  FD  method  has  its  worst  behavior  with  UD,  used  in 
[39],  as  shown  in  Fig.  2a. 

7.3.  Frequency  domain  tests 

Results  discussed  so  far  give  a  general  outlook  about  the  accu¬ 
racy  of  different  discretization  methods  along  with  some  grids. 
However,  the  test  function  used  is  very  specific  and  similar  results 
cannot  be  ensured  in  general  for  different  input  profiles.  Thus,  dis¬ 
cretization  methods  were  tested  also  in  the  frequency  domain.  In 
this  case,  just  the  best  and  worst  approximations  distinguished 
with  the  time  domain  tests  were  evaluated.  Approximations 
obtained  with  PDQand  FDQalong  to  PD,  the  worst  scenario  of  these 
methods,  were  also  tested  in  the  frequency  domain  despite  they 
are  unfeasible  because  error  calculated  in  time  domain  worsens  as 
order  q  increases. 

Frequency  response  of  magnitude  and  phase  of  worst  approx¬ 
imations  are  shown  in  Fig.  4a  and  b,  whereas  responses  of  best 
approximations  are  shown  in  Fig.  4c  and  d.  In  both  pairs  of  fig¬ 
ures  the  base  response  is  that  of  the  analytic  TF  (42).  Responses 
of  approximations  (50)  and  (51 ),  which  can  also  be  brought  to  the 
general  form  (14)  and  (15),  are  also  included  in  figures.  As  similar 


results  were  observed  for  higher  and  lower  orders,  semidiscretiza¬ 
tion  approximations  of  order  q  =  14  were  arbitrarily  selected  as 
examples  for  frequency  response  tests  (14  is  the  order  for  approx¬ 
imations  of  the  positive  electrode  subsystem  when  truncation  is 
based  on  bandwidth  selection).  Also,  in  Fig.  4a-d,  the  vertical  bar 
represents  the  normalized  module  of  the  faster  analytic  pole  pq  of 
approximations,  mentioned  in  Section  5  and,  according  to  (40),  in 
this  case  is  |pi4//<j  =  X\4  =  1796.7.  For  more  generality,  the  nor¬ 
malized  frequency  oo  =  oo/K  is  defined. 


0.35  - 
0.3  - 
0.25  - 
0.2  - 
0.15  - 
0.1 

0.05  - 
0  - 


[ 

] 

c 

rtf- 

4  5  6  7 

Index  of  residues  m 


Fig.  5.  Normalized  residues  associated  to  the  pulsating  input  signal  ft. 
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a  Cell  voltage  using  low-order  polynomial-based  approximations 


b  Negative  electrode  state-of-charge  using  low-order  polynomial-based 
approximations 


C  Cell  voltage  using  semidiscretization  approximations  and  truncation  d  Negative  electrode  state-of-charge  using  semidiscretization  approxi- 


by  bandwidth  selection 


mations  and  truncation  by  bandwidth  selection 


e  Cell  voltage  using  semidiscretization  approximations  and  truncation  26  f  Negative  electrode  state-of-charge  using  semidiscretization  approxi- 
by  residue  analysis  mations  and  truncation  by  residue  analysis 

Fig.  6.  Results  of  test  applying  the  HPPC  current  profile. 


Because  of  their  low  order,  approximations  (50)  and  (51) 
have  poor  but  reasonable  responses.  Particularly,  the  magnitude 
response  of  (51),  used  in  [33],  is  almost  as  good  as  that  of  the 
FD  14th-order  approximations  when  using  UD.  However,  the 
phase  responses  of  (50)  and  (51)  shows  the  limitations  of  such 
approximations,  being  valid  for  low  frequencies  in  the  interval  of 
boe  [0,  10]. 

The  worst  approximations  (see  Fig.  4a  and  b)  of  PDQ  and 
FDQ  showed  to  have  similar  magnitude  responses,  better  to  that 


attached  to  FD  and  FE  despite  their  growing  error  obtained  in  the 
time  tests.  PDQand  FDQ  approximations  sustain  a  good  agreement 
with  the  magnitude  base  response  until  half  decade  after  bo  =  ^14 
and  FD  moves  away  almost  one  decade  before  reaching  X\4.  On  the 
other  hand,  the  phase  responses  of  FD  diverge  too  early  from  the 
base  response  and  that  of  PDQ  and  FDQ  keep  their  agreement  to 
the  base  phase  response  until  X\4.  Paradoxically,  RLP  distribution 
gives  the  worst  performance  to  FE  in  the  time  domain,  but  causes 
its  best  frequency  response. 
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On  the  other  hand,  in  Fig.  4c  and  d,  it  is  shown  that  the  response 
of  semidiscretization  approximations  improves  when  using  the 
grids  which  also  drove  to  the  best  results  in  the  time  domain.  PDQ 
and  FDQ  magnitude  response  is  close  to  the  base  response  up  to 
one  decade  after  A24  and  half  decade  for  phase  response.  It  should 
be  highlighted  that  behavior  of  FD  improves  when  using  RLP  to  the 
extent  that  such  approximations  is  almost  as  good  as  PDQ  in  both 
time  and  frequency  domains.  Flowever,  for  FE  the  result  is  opposite 
again;  the  frequency  response  worsens  with  PD. 

7.4.  Test  with  the  HPPC  signal 

The  response  to  the  FIPPC  signal  of  the  single-particle 
electrochemical  model  is  presented.  The  voltage  evolution  of 
approximations  against  the  response  profile  reported  in  [35]  for 
58.3%  of  battery  SOC,  acquired  from  experimental  tests  in  the  same 
work,  is  shown.  For  tests,  initial  SOC  was  slightly  adjusted  to  57%  as 
well  as  Rf  from  20  to  23  £2  cm2 .  The  pair  of  fast  peaks  at  52  s  was  not 
induced  because  of  the  lack  of  information  about  them.  Neverthe¬ 
less,  due  to  the  sluggish  nature  of  the  system,  their  influence  to  the 
cell  response  is  not  substantial.  The  evolution  of  negative  electrode 
state-of-charge  0n  is  also  shown  as,  for  this  electrode,  results  from 
truncation  criteria  are  more  drastic.  In  this  case,  the  base  solution 
was  solved  numerically  using  the  pdepe  routine  introduced  above. 
Fig.  6a  and  b  correspond  to  cell  response  when  approximations 
(50)  and  (51)  are  applied.  It  is  reported  in  the  literature  that  those 
approximations  work  well  for  constant  long  time  and  low  input  cur¬ 
rent  regimes  (see  for  example  [20,28])  but,  as  shown,  fail  dramati¬ 
cally  for  pulsating  profiles  due  to  their  poor  bandwidth  and  thus  is 
not  advisable  to  use  this  approach  for  complicated  input  signals. 

For  spatial  semidiscretization  approximations,  the  most  typical 
grid  was  considered  for  each  case.  UD  is  used  for  FD  as  in  [39]. 
PD  is  used  for  FE,  which  is  similar  to  the  grid  used  in  [35].  Finally, 
CGL  is  used  along  with  PDQ  and  FDQ  as  recommended  in  [44].  The 
response  of  approximations  is  quite  good  when  model  order  is  cho¬ 
sen  by  bandwidth  selection.  Voltage  evolution  of  approximations 
is  shown  in  Fig.  6c.  Approximations  related  to  FE,  PDQ  and  FDQ 
methods  follow  very  closely  the  experimental  data  taken  from  [35]. 
In  particular,  PDQ  and  FDQ  are  better  than  FE  and  seem  to  have 
both  the  same  good  accuracy.  Flowever,  FD  used  along  with  UD  as 
reported  in  the  literature,  produces  a  noticeable  error  because  it  is 
the  worst  grid  choice  for  that  method.  Nevertheless,  as  discussed 
before,  the  error  of  FD  can  be  improved  whether  RLP  and  CGL  grids 
are  used  instead  of  UD.  In  Fig.  6d,  error  of  approximations  is  more 
evident  than  in  Fig.  6c  because  in  the  latter  case  0n  is  affected  by 
function  (A.2)  of  the  negative  electrode  open  circuit  potential. 

The  response  of  the  single-particle  model  when  the  order  of 
approximations  is  chosen  through  the  residue  analysis  is  shown 
in  Fig.  6e  and  f.  Voltage  responses  for  PDQ  and  FDQ  seem  not  to 
undergo  important  loss  of  quality  despite  the  model  order  is  less 
than  half  of  that  obtained  by  bandwidth  selection.  Indeed,  in  this 
situation,  these  methods  offer  a  better  voltage  response  than  FD 
when  order  is  defined  by  bandwidth  selection.  The  response  of 
FE  approximation  slightly  worsens  with  the  order  reduction  and 
error  of  FD  increases  in  a  great  extent.  This  is  clearer  again  in  the 
response  of  negative  electrode  state-of-charge.  Whereas  PDQ  and 
FDQ  approximations  could  be  still  used  for  state-of-charge  estima¬ 
tion,  FD  and  FE  approximations  are  not  useful  when  order  is  chosen 
by  the  second  criterion. 

8.  Conclusions 

In  this  paper  the  state  equation  of  the  single-particle  elec¬ 
trochemical  model  was  studied  with  the  goal  of  finding  reduced 
order  approximations  that  represent  adequately  the  battery 
main  dynamics  for  some  given  operation  conditions.  Spatial 


semidiscretization  approach  was  reintroduced  for  order  reduc¬ 
tion  and  three  discretization  methods,  combined  with  four  grids, 
were  applied  in  a  classical  fashion:  finite  difference,  finite  element 
and  differential  quadrature.  Because  of  their  recurrent  appear¬ 
ance  in  the  literature,  low  order  dynamical  approximations  based 
on  second  and  fourth  order  polynomials  were  also  evaluated  and 
compared  against  semidiscretization  approximations.  Addition¬ 
ally,  two  order  selection  or  truncation  criteria,  suggested  but  not 
explored  in  the  literature,  were  discussed  and  contrasted.  The  first, 
more  conservative,  consists  on  selecting  the  approximation  band¬ 
width  such  that  it  covers  at  least  the  major  part  of  the  energy 
spectrum  of  some  proposed  signal  test.  The  second  criterion  is 
based  in  a  residue  analysis  of  the  convolution  between  the  state 
equation  and  the  test  signal  and  drove  to  more  practical  results. 

Approximations  were  first  evaluated  in  the  time  domain  to 
have  a  preliminary  surmise  about  which  had  the  best  trade  off 
between  low  order  and  accuracy.  In  this  case,  the  test  signal  con¬ 
sisted  in  a  rectangular  pulse  followed  by  a  relaxation  time,  whose 
duration  was  such  that  allowed  to  observe  the  whole  transient 
response  under  forced  as  well  as  under  free  excitation.  The  fre¬ 
quency  response  of  approximations  was  also  compared  against 
that  of  the  analytic  transfer  function  related  to  the  state  equation 
in  order  to  confirm  the  results  obtained  in  the  time  domain  tests. 
Finally,  representative  cases  of  all  considered  approximations  were 
applied  to  the  single-particle  model  and  their  response  to  the  FIPPC 
test  signal  was  compared  against  experimental  data  extracted  from 
the  literature.  Best  results  were  obtained  with  discretizations  based 
on  polynomial  differential  quadrature  whereas  it  was  also  observed 
that  approximations  based  on  low  order  polynomials  have  a  very 
poor  performance  for  pulsating  input  currents. 
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Table  B.l 

Values  of  single-particle  model  parameters  and  constants. 


Descripion 

Symbol 

Value 

Solid  phase  Li+  diffusivity 

Dn 

2  x  10“12  cm2  s_1 

DP 

3.7  x  10~12  cm2  s”1 

Representative  particle  radius 

Rn 

1  x  10“4  cm 

RP 

1  x  10“4  cm 

Specific  contact  area a 

an 

17400  cm-1 

ap 

15  000  cm-1 

Element  thickness 

Sn 

50  x  10-4  cm 

$sep 

25.4  x  10“4  cm 

Sp 

36.4  x  10“4  cm 

Sectional  cell  area 

A 

10452  cm2 

Film  resistance  at  interphase 

Rf 

20  £2  cm2 

Effective  solution  conductivity  a 

Kn 

10.87  xlO“3S cm-1 

Ksep 

20.08  xlO“3S cm-1 

Kp 

10.77  xlO“3S cm-1 

Maximum  Li+  concentration 

CnmaX 

16.1  x  10“3  mol  cm-3 

rmax 

P 

23.9  x  10“3  mol  cm-3 

Exchange  current  density 

lo,n 

3.6  x  10“3  A  cm-2 

io,p 

2.6  x  10-3  A  cm-2 

Stoichiometry  at  0%  SOC 

6>° 

0.126 

0°P 

0.936 

Stoichiometry  at  100%  SOC 

<9™ 

0.676 

or 

0.442 

Subscripts 

n 

Negative  electrode 

sep 

Separator 

P 

Positive  electrode 

Faraday’s  constant 

T 

96  487  C  mol-1 

Ideal  gas  constant 

n 

8.3143  J  mol”1  K-1 

Absolute  temperature 

T 

293.15  K  (20  °C) 

a  Calculated  from  data  in  [29]. 


A.  Romero-Becerril,  L.  Alvar ez-Icaza  /  Journal  of  Power  Sources  196  (201 1 )  10267-10279 


10279 


Appendix  A.  Open  circuit  potentials 

For  the  positive  and  negative  electrodes,  respectively: 
l/p(0p)  =  85.6810®  -375.700® +  613.8904-555.6503  +281.0602 
—  76.6480p  +  13.1983  -  0.30987 exp  (5.6570”®)  (A.l) 


U„(0n)  =  8.00229  +  5.06470,,  -  12.578O0y2  -  8.6322  x  1O"40-1 
+  2.1765  x  1O“502/2  -0.46016exp[15. 0(0.06 -0„)] 

-  0.55364 exp[-2.4326(0„  -0.92)]  (A.2) 


Appendix  B.  List  of  parameters 

See  Appendix  Table  B.l 
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